#########
# Setup #
#########

# Instantiate
	library(foreign)
	rueda = read.csv("Duplication dataset.csv")
	attach(rueda)

# kill bogus columns
	rueda$X<- NULL
	rueda$X.1<- NULL

# rename minwage
	minwag = rueda$min

# Libraries

	library("lmtest")
	library("plm")
	library("apsrtable")
	library("xtable")


#############
#  Table 1 	#
#############

m5010 = c()
diff=c()
for (i in 1:max(count)){
 ## means for each country
 m5010[i] = mean(ia5010[count==i],na.rm=T)

 ## % change
 firstyrineq = ia5010[count==i & year==min(year[ia5010!="NA" & count==i], na.rm=T)]
 secyrineq = ia5010[count==i & year==max(year[ia5010!="NA" & count==i], na.rm=T)]
 diff[i] = (secyrineq - firstyrineq)*100 / firstyrineq
}

CountryName = unique(country)
Mean = round(m5010,2)
Delta = round(diff,2)
averages = cbind( round(mean(m5010),2), round(mean(diff),2) )
sds = cbind( round(sd(m5010),2),round(sd(diff),2) )

r_estimates <- data.frame(CountryName,Mean,Delta)
r_average <- data.frame("Average:",averages)
r_sd <- data.frame("Standard Deviation:",sds)
colnames(r_average) <- colnames(r_estimates)
colnames(r_sd) <- colnames(r_estimates)

results <- rbind(r_estimates,r_average,r_sd)
colnames(results) <- c("Country", "Mean", "Delta")

library(xtable)
newobject<-xtable(results)
print(newobject, floating = FALSE)


##############
## Table 3 	 #
##############

varsofint = rbind(govpart,hkcorp,govem,generos,minwag)
vmeans = rowMeans(varsofint,na.rm=T)

mn=c()
mx=c()
std=c()
for(i in 1:5){
 std[i] = sd(varsofint[i,],na.rm=T)
 mn[i] = min(varsofint[i,],na.rm=T)
 mx[i] = max(varsofint[i,],na.rm=T)
}

results <- cbind(vmeans,std,mn,mx)
colnames(results) <- c("Mean","Standard Deviation", "Minimum", "Maximum")
rownames(results) <- c("Cabinet Partisanship","Corporatism", "Government Employment", "Welfare State Generosity","Minimum Wages")
newobject<-xtable(results)
print(newobject, floating = FALSE)


##########################
## Appendix 2/Figure7 	 #
##########################


fit=lm(govem~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90+factor(country),data=rueda)
dvemp2 <- plm(govem~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")
test <- coeftest(dvemp2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
append = c(1:4)
test <- rbind(append,test)
fit$se=test[,3]
fit$coef=test[,1]

fit2=lm(generos~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90+factor(country),data=rueda)
dvwelf2 = plm(generos~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")
test <- coeftest(dvwelf2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
append = c(1:4)
test <- rbind(append,test)
fit2$se=test[,3]
fit2$coef=test[,1]

fit3=lm(minwag~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90+factor(country),data=rueda)
dvmin2 = plm(minwag~govpart+hkcorp+govpart:hkcorp+open+finop+dccggx+unemp+gdpgr+deca70+deca80+deca85+deca90,data=rueda,model="within")
test <- coeftest(dvmin2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
append = c(1:4)
test <- rbind(append,test)
fit3$se=test[,3]
fit3$coef=test[,1]

apsrtable(fit, fit2, fit3, digits=3,omitcoef=list("(Intercept)"),se="robust",stars="default") 


dvemp3 = coeftest(dvemp2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
dvwelf3 = coeftest(dvwelf2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
dvmin3 = coeftest(dvmin2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))

# Marginal effects 

  lowcorp = .15
	highcorp = .9


# Partisanship --> Government Employment

	# Low corp:
	dvemp2$coeff[1] + (lowcorp*dvemp2$coeff[12])
		# = 0.03250646

      (dvemp2$coeff[1] + (lowcorp*dvemp2$coeff[12])) / sqrt( dvemp3["govpart",2]^2 +(lowcorp^2)*dvemp3["govpart:hkcorp",2]^2 + 2*lowcorp*( vcovBK(dvemp2, diagonal=T,cluster="time")[1,12]) ) 
            # tstat = 2.809629 

	# High corp  
	dvemp2$coeff[1] + (highcorp*dvemp2$coeff[12])
		#=  -0.0384865
      (dvemp2$coeff[1] + (highcorp*dvemp2$coeff[12])) / sqrt( dvemp3["govpart",2]^2 +(highcorp^2)*dvemp3["govpart:hkcorp",2]^2 + 2*highcorp*( vcovBK(dvemp2, diagonal=T,cluster="time")[1,12]) ) 
            # tstat = -2.667466 

# Partisanship --> Welfare Generosity

	# Low corp:
	dvwelf2$coeff[1] + (lowcorp*dvwelf2$coeff[12])
		 ## -0.02006564 
      (dvwelf2$coeff[1] + (lowcorp*dvwelf2$coeff[12])) / sqrt( dvwelf3["govpart",2]^2 +(lowcorp^2)*dvwelf3["govpart:hkcorp",2]^2 + 2*lowcorp*( vcovBK(dvwelf2, diagonal=T,cluster="time")[1,12]) ) 
            # -1.268140 

	# High corp  
	dvwelf2$coeff[1] + (highcorp*dvwelf2$coeff[12])
	 ## -0.06432758 
     (dvwelf2$coeff[1] + (highcorp*dvwelf2$coeff[12])) / sqrt( dvwelf3["govpart",2]^2 +(highcorp^2)*dvwelf3["govpart:hkcorp",2]^2 + 2*highcorp*( vcovBK(dvwelf2, diagonal=T,cluster="time")[1,12]) ) 
            # -2.504921 

# Partisanship --> Minimum Wage

	# Low corp:
	dvmin2$coeff[1] + (lowcorp*dvmin2$coeff[12])
		 ## 0.001047811
      (dvmin2$coeff[1] + (lowcorp*dvmin2$coeff[12])) / sqrt( dvmin3["govpart",2]^2 +(lowcorp^2)*dvmin3["govpart:hkcorp",2]^2 + 2*lowcorp*( vcovBK(dvmin2, diagonal=T,cluster="time")[1,12]) ) 
            # 2.858847 

	# High corp  
	dvmin2$coeff[1] + (highcorp*dvmin2$coeff[12])
	## -0.000487458 
     (dvmin2$coeff[1] + (highcorp*dvmin2$coeff[12])) / sqrt( dvmin3["govpart",2]^2 +(highcorp^2)*dvmin3["govpart:hkcorp",2]^2 + 2*highcorp*( vcovBK(dvmin2, diagonal=T,cluster="time")[1,12]) ) 
            # -2.231917


fit=lm(ia5010~govem+generos+minwag+hkcorp+govem:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90+factor(country),data=rueda)
ivemp2 = plm(ia5010~govem+generos+minwag+hkcorp+govem:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")
test <- coeftest(ivemp2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
append = c(1:4)
test <- rbind(append,test)
fit$se=test[,3]
fit$coef=test[,1]

fit2=lm(ia5010~govem+generos+minwag+hkcorp+generos:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90+factor(country),data=rueda)
ivwelf2 = plm(ia5010~govem+generos+minwag+hkcorp+generos:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")
test <- coeftest(ivwelf2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
append = c(1:4)
test <- rbind(append,test)
fit2$se=test[,3]
fit2$coef=test[,1]

fit3=lm(ia5010~govem+generos+minwag+hkcorp+minwag:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90+factor(country),data=rueda)
ivmin2 = plm(ia5010~govem+generos+minwag+hkcorp+minwag:hkcorp+unemp+ldc+femla+priv+deca70+deca80+deca85+deca90,data=rueda,model="within")
test <- coeftest(ivmin2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
append = c(1:4)
test <- rbind(append,test)
fit3$se=test[,3]
fit3$coef=test[,1]

apsrtable(fit, fit2, fit3, digits=3, omitcoef=list("(Intercept)"),se="robust",stars="default") 


ivemp3 = coeftest(ivemp2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
ivwelf3 = coeftest(ivwelf2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
ivmin3 = coeftest(ivmin2, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))

### Marginal effects

# Government Employment -> Inequality

	# Low corp:
	ivemp2$coeff[1] + (lowcorp*ivemp2$coeff[13])
		# =  -0.01119048
	(ivemp2$coeff[1] + (lowcorp*ivemp2$coeff[13])) / sqrt( ivemp3["govem",2]^2 +(lowcorp^2)*ivemp3["govem:hkcorp",2]^2 + 2*lowcorp*( vcovBK(ivemp2, diagonal=T,cluster="time")[1,13]) ) 
     		# - 2.824377

	# High corp  
	ivemp2$coeff[1] + (highcorp*ivemp2$coeff[13])
		#=  -0.008057232 
	(ivemp2$coeff[1] + (highcorp*ivemp2$coeff[13])) / sqrt( ivemp3["govem",2]^2 +(highcorp^2)*ivemp3["govem:hkcorp",2]^2 + 2*highcorp*( vcovBK(ivemp2, diagonal=T,cluster="time")[1,13]) ) 
		# 2.596

#  Welfare Generosity -> Inequality

	# Low corp:
	ivwelf2$coeff[2] + (lowcorp*ivwelf2$coeff[13])
		# = 0.002068815 
 	(ivwelf2$coeff[2] + (lowcorp*ivwelf2$coeff[13])) / sqrt( ivwelf3["generos",2]^2 +(lowcorp^2)*ivwelf3["generos:hkcorp",2]^2 + 2*lowcorp*( vcovBK(ivwelf2, diagonal=T,cluster="time")[2,13]) ) 
		# .924

	# High corp  
	ivwelf2$coeff[2] + (highcorp*ivwelf2$coeff[13])
		#=  -0.0014279 
      (ivwelf2$coeff[2] + (highcorp*ivwelf2$coeff[13])) / sqrt( ivwelf3["generos",2]^2 +(highcorp^2)*ivwelf3["generos:hkcorp",2]^2 + 2*highcorp*( vcovBK(ivwelf2, diagonal=T,cluster="time")[2,13]) ) 
		# -1.188

# Minimum Wage -> Inequality

	# Low corp:
	ivmin2$coeff[3] + (lowcorp*ivmin2$coeff[13])
		# = -0.2519744
	(ivmin2$coeff[3] + (lowcorp*ivmin2$coeff[13])) / sqrt( ivmin3["minwag",2]^2 +(lowcorp^2)*ivmin3["minwag:hkcorp",2]^2 + 2*lowcorp*( vcovBK(ivmin2, diagonal=T,cluster="time")[3,13]) ) 
		# -2.4179

	# High corp  
	ivmin2$coeff[3] + (highcorp*ivmin2$coeff[13])
		#=  0.2227287
 	(ivmin2$coeff[3] + (highcorp*ivmin2$coeff[13])) / sqrt( ivmin3["minwag",2]^2 +(highcorp^2)*ivmin3["minwag:hkcorp",2]^2 + 2*highcorp*( vcovBK(ivmin2, diagonal=T,cluster="time")[3,13]) ) 
		# 1.014


################################
## Figure8 ##
################################



####### Policies as DV

controls = cbind(open,finop,dccggx,unemp,gdpgr)
policydv = cbind(govem,generos,minwag)

vlowcorp=matrix(NA,3,5) 
vhighcorp=matrix(NA,3,5)
tstatlowcorp = matrix(NA,3,5)
tstathighcorp = matrix(NA,3,5)

# Outer loop - dependent variables
for(j in 1:3){ 

# Inner loop - IVs
	for(i in 1:5){
  		exreg = plm(policydv[,j]~govpart+hkcorp+govpart:hkcorp+controls[,-i]+deca70+deca80+deca85+deca90,data=rueda,model="within")
 	      exreg1 = coeftest(exreg, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
  		## marginal effects
 		 	vlowcorp[j,i] = exreg1["govpart",1] + lowcorp*exreg1["govpart:hkcorp",1]
  		 	vhighcorp[j,i]= exreg1["govpart",1] + highcorp*exreg1["govpart:hkcorp",1]
 		## t-stats for marginal effects
               tstatlowcorp[j,i] = (exreg1["govpart",1] + lowcorp*exreg1["govpart:hkcorp",1])/ sqrt( exreg1["govpart",2]^2 +(lowcorp^2)*exreg1["govpart:hkcorp",2]^2 + 2*lowcorp*(vcovBK(exreg, diagonal=T,cluster="time")[1,11]) )
 	         tstathighcorp[j,i] = (exreg1["govpart",1] + highcorp*exreg1["govpart:hkcorp",1])/ sqrt( exreg1["govpart",2]^2 +(highcorp^2)*exreg1["govpart:hkcorp",2]^2 + 2*highcorp*(vcovBK(exreg, diagonal=T,cluster="time")[1,11]) )
          }
}

results=matrix(NA,6,5)
for (i in 1:5){
	results[1,i] = round(vlowcorp[1,i],4)
	results[2,i] = round(tstatlowcorp[1,i],2)
	results[3,i] = round(vlowcorp[2,i],4)
	results[4,i] = round(tstatlowcorp[2,i],2)
	results[5,i] = round(vlowcorp[3,i],4)
	results[6,i] = round(tstatlowcorp[3,i],2)
}
## this results matrix gives all estimates and t-stats at low corp
## the bounds can be read off it
results

results=matrix(NA,6,5)
for (i in 1:5){
	results[1,i] = round(vhighcorp[1,i],4)
	results[2,i] = round(tstathighcorp[1,i],2)
	results[3,i] = round(vhighcorp[2,i],4)
	results[4,i] = round(tstathighcorp[2,i],2)
	results[5,i] = round(vhighcorp[3,i],4)
	results[6,i] = round(tstathighcorp[3,i],2)
}
## equivalent for high corp
results


####### Policies as IV

controls2 = cbind(unemp,ldc,femla,priv)
policyiv = cbind(govem,generos,minwag)

vlowcorp=matrix(NA,3,4) 
vhighcorp=matrix(NA,3,4)
tstatlowcorp=matrix(NA,3,4) 
tstathighcorp=matrix(NA,3,4)

for(j in 1:3){

 for(i in 1:4){
		exreg = plm(ia5010~govem+generos+minwag+hkcorp+policyiv[,j]:hkcorp+controls2[,-i]+deca70+deca80+deca85+deca90,data=rueda,model="within")
 	      exreg1 = coeftest(exreg, vcov=function(x) vcovBK(x, diagonal=T,cluster="time"))
		## marginal effects
            vlowcorp[j,i] = exreg$coeff[j] + lowcorp*exreg$coeff[length(exreg$coeff)]
  		vhighcorp[j,i]= exreg$coeff[j] + highcorp*exreg$coeff[length(exreg$coeff)]
            ## t-stats for marginal effects
            tstatlowcorp[j,i] = vlowcorp[j,i] / sqrt( exreg1[j,2]^2 +(lowcorp^2)*exreg1[nrow(exreg1),2]^2 + 2*lowcorp*(vcovBK(exreg, diagonal=T,cluster="time")[j,length(exreg$coeff)]) )
 	      tstathighcorp[j,i] = vhighcorp[j,i] / sqrt( exreg1[j,2]^2 +(highcorp^2)*exreg1[nrow(exreg1),2]^2 + 2*highcorp*(vcovBK(exreg, diagonal=T,cluster="time")[j,length(exreg$coeff)]) )
   }
      
}


results=matrix(NA,6,4)
for (i in 1:4){
	results[1,i] = round(vlowcorp[1,i],4)
	results[2,i] = round(tstatlowcorp[1,i],2)
	results[3,i] = round(vlowcorp[2,i],4)
	results[4,i] = round(tstatlowcorp[2,i],2)
	results[5,i] = round(vlowcorp[3,i],4)
	results[6,i] = round(tstatlowcorp[3,i],2)
}
## this results matrix gives all estimates and t-stats at low corp
## the bounds can be read off it
results

results=matrix(NA,6,4)
for (i in 1:4){
	results[1,i] = round(vhighcorp[1,i],4)
	results[2,i] = round(tstathighcorp[1,i],2)
	results[3,i] = round(vhighcorp[2,i],4)
	results[4,i] = round(tstathighcorp[2,i],2)
	results[5,i] = round(vhighcorp[3,i],4)
	results[6,i] = round(tstathighcorp[3,i],2)
}
## equivalent for high corp
results
